#!/usr/bin/perl -w
use strict;


package pack_ascore_TMT;


sub ascore{
	my ($peptide,$parent_mass,$z,$masses,$intensities,$precision,$limit_mz_sup,$limit_mz_inf)=@_;
	#print "peptide standard notation: $peptide\n parent: $parent_mass \n charge: $z \n $masses \n $intensities \n";
	#print $peptide,"\n";
	my %result;
	#El numero de sitios de fosforilacion ha de ser mayor que el de 
	#fosforilaciones para que el Ascore tenga sentido. Si el Ascore no tiene sentido, lo devuelvo como 1000
		my $numero_sitios_fosforilacion=0;
		my $numero_sitios_deshidratacion=0;
		my $numero_sitios_modificacion=0;
		my $numero_fosforilaciones=0;
		my $numero_deshidrataciones=0;
		my $numero_modificaciones=0;
		
		my $peptido_conteo_modificaciones=$peptide;
		while($peptido_conteo_modificaciones=~s/(\(21\)|\(21\,|21\))/x/){
			$numero_fosforilaciones++;
		}
		while($peptido_conteo_modificaciones=~s/(\(23\)|\(23\,|23\))/x/){
			$numero_deshidrataciones++;
		}
		while($peptido_conteo_modificaciones=~s/^Y/x/){
			$numero_sitios_fosforilacion++;
			$numero_sitios_modificacion++;
		}
		while($peptido_conteo_modificaciones=~s/^[ST]/x/){
			$numero_sitios_fosforilacion++;
			$numero_sitios_deshidratacion++;
			$numero_sitios_modificacion++;
		}
		while($peptido_conteo_modificaciones=~s/Y(\w)/x$1/){
			$numero_sitios_fosforilacion++;
			$numero_sitios_modificacion++;
		}
		while($peptido_conteo_modificaciones=~s/[ST](\w)/x$1/){
			$numero_sitios_fosforilacion++;
			$numero_sitios_deshidratacion++;
			$numero_sitios_modificacion++;
		}

	#Si el numero de sitios de modificacion es el mismo que el de modificaciones
	$numero_modificaciones=$numero_fosforilaciones+$numero_deshidrataciones;
	if($numero_sitios_modificacion==$numero_modificaciones){
		$result{'phosphorylations'}=$numero_fosforilaciones;
		$result{'dehydratations'}=$numero_deshidrataciones;
		$result{'peptide'}=$peptide;
		if($numero_modificaciones==1){
			$result{'ascore'}=[1000];
		}
		elsif($numero_modificaciones==2){
			$result{'ascore'}=[1000,999];
		}
		elsif($numero_modificaciones==3){
			$result{'ascore'}=[1000,999,998];
		}
		return %result;
	}
	#Caso especial: numero de deshidrataciones es igual al numero de posibles sitios de deshidratacion
		if($numero_deshidrataciones>0){
			if($numero_sitios_deshidratacion==$numero_deshidrataciones){
				if($numero_modificaciones==1){
					$result{'phosphorylations'}=$numero_fosforilaciones;
					$result{'dehydratations'}=$numero_deshidrataciones;
					$result{'peptide'}=$peptide;
					$result{'ascore'}=[1000];
					return %result;
				}
				else{
					#Caso no resuelto: num.DHs == num.sitios.DH Y num_mods>=2
					$result{'phosphorylations'}=$numero_fosforilaciones;
					$result{'dehydratations'}=$numero_deshidrataciones;
					$result{'peptide'}=$peptide;
					if($numero_modificaciones==2){
						$result{'ascore'}=[0,0];
					}
					elsif($numero_modificaciones==3){
						$result{'ascore'}=[0,0,0];
					}
					print "CASO NO RESUELTO \n";
				}
				
			}
		}
	if($numero_modificaciones==0){
		$result{'phosphorylations'}=$numero_fosforilaciones;
		$result{'dehydratations'}=$numero_deshidrataciones;
		$result{'peptide'}=$peptide;
		return %result;
	}
	elsif($numero_modificaciones==1){
		%result=ascore_1_fosforilacion($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
		return %result;
	}
	elsif($numero_modificaciones==2){
		%result=ascore_2_fosforilaciones($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
		return %result;
	}
	elsif($numero_modificaciones==3){
		%result=ascore_3_fosforilaciones($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
		return %result;
	}
	else{
		$result{'phosphorylations'}=$numero_fosforilaciones;
		$result{'dehydratations'}=$numero_deshidrataciones;
		$result{'peptide'}=$peptide;
		return %result;
	}
	
	

}

sub ascore_1_fosforilacion{
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my $masses=$_[2];			#referencia a array de masas
	my $intensities=$_[3];	#referencia a array intensidades
	my $z=$_[4];				#carga
	my $limit_mz_sup=$_[5];			
	my $limit_mz_inf=$_[6];
	#print "Ascore 1 fosfo\n";
	my @masas=split /\|/,$masses;
	my @intensidades=split /\|/,$intensities;
	

	#print $peptido,"\n";#exit;
	my @peptidos_posibles=peptide_collection($peptido);
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z,$limit_mz_sup,$limit_mz_inf)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{'score'}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{'score_ponderado'}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{'score'}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}
	
	#Ordeno peptidos segun Peptide score
	my @peptidos_con_sumas_scores_ordenadas=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_sumas_scores_ordenadas[0];
	my $p2=$peptidos_con_sumas_scores_ordenadas[1];
	my $p3=$peptidos_con_sumas_scores_ordenadas[2];
	
	#Convierto p1 y p2 a peptidos en formato [styuv]
	my $p1_s=conversion_styuv($p1);
	my $p2_s=conversion_styuv($p2);

	#print $p1_s,"\n",$p2_s,"\n";
	#if(($peptide_scores{$p2})and($peptide_scores{$p1}==$peptide_scores{$p2})){
		#print "Peptide scores iguales para peptidos 1 y 2\n";
		#return ($peptido,$peptide_scores{$p1},0,0);
	#}
	#testeo si, en el caso de haber un tercer candidato, su peptide score es tan bueno como el segundo;
	if (defined $p3){ 
		my $p3_s=conversion_styuv($p3);			
		if((defined $peptide_scores{$p3})and($peptide_scores{$p2}==$peptide_scores{$p3})){
			#Si p2 tiene el mismo Peptide Score que p3, debo elegor uno para usar en el ascore 
			#print "pep2 ($p2) y pep 3 ($p3) son iguales ($peptide_scores{$p2}), respecto a $p1\n";
			###Determino cual de los peptidos p2 y p3 tiene la fosforilacion mas proxima a p1
			my @p_trozos_peptido_1= split //,$p1_s;
			my @p_trozos_peptido_2= split //,$p2_s;
			my @p_trozos_peptido_3= split //,$p3_s;
			my ($punto1,$punto2,$punto3);
			for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){if($p_trozos_peptido_1[$i]=~/[styuv]/){$punto1=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){if($p_trozos_peptido_2[$i]=~/[styuv]/){$punto2=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_3+1;$i++){if($p_trozos_peptido_3[$i]=~/[styuv]/){$punto3=1+$i;}}
			my $p2_usar; #uso el peptido con menor diferencia a punto1
			$punto2=abs($punto1-$punto2);
			$punto3=abs($punto1-$punto3);
			if($punto2<$punto3){
				$p2_usar=$p2;
			}
			else{
				$p2_usar=$p3;
			}
			$p2=$p2_usar;
			#print "Peptido usado: ".$p2_usar."\n";
			
		}
	}
	
	my @peptidos_ascore;
	push @peptidos_ascore,$p1;
	push @peptidos_ascore,$p2;
	

	#Obtengo nivel maxima diferencia entre Peptide Scores
	my $nivel_maxima_diferencia=0;
	my $difss=0;
	for(my $i=1;$i<11;$i++){
		my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$p2}{$i}{score};
		if($dif>$difss){
			$difss=$dif;
			$nivel_maxima_diferencia=$i;
		}
	}
	if($nivel_maxima_diferencia==0){
		$nivel_maxima_diferencia=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
	}
	
	my @p_trozos_peptido_1= split //,$p1_s;
	my @p_trozos_peptido_2= split //,$p2_s;
	my ($punto1,$punto2);
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			$punto1=1+$i;
		}
	}
	for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){
		if($p_trozos_peptido_2[$i]=~/[styuv]/){
			$punto2=1+$i;
		}
	}
	if($punto1>$punto2){
		($punto1,$punto2)=($punto2,$punto1);
	}
	
	my %fragmentos_p1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p1=%{tratamiento_fragmentos($p1_s,\%fragmentos_p1,$punto1,$punto2)};
	my %fragmentos_p2=%{prever_fragmentos($p2,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p2=%{tratamiento_fragmentos($p2_s,\%fragmentos_p2,$punto1,$punto2)};

	
	my %datos_ascore;
	$datos_ascore{$p1}=matches_ascore(\%fragmentos_p1,\@masas,\@intensidades,$precision,$p1);
	$datos_ascore{$p2}=matches_ascore(\%fragmentos_p2,\@masas,\@intensidades,$precision,$p2);
			
	my $ascore;
	#print "p1-> $datos_ascore{$p1}{$nivel_maxima_diferencia}{score}\n";
	#print "p2-> $datos_ascore{$p2}{$nivel_maxima_diferencia}{score}\n";
	
	my $matches1=1+$#{$datos_ascore{$p1}{$nivel_maxima_diferencia}{tags_match_intervalo}};
	my $matches2=1+$#{$datos_ascore{$p2}{$nivel_maxima_diferencia}{tags_match_intervalo}};
	my @array1=keys(%fragmentos_p1);
	my $size=1+$#array1;
	my $proporcion=$matches1."//".$size." (".$matches2.")";
	
	#if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})and(not defined $datos_ascore{$p2}{$nivel_maxima_diferencia}{score})){
		#print "No definido score p1 ni score p2\n";
		#return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
	#}
	#if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})){
		#print "No definido score p1\n";
		#return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
	#}

	$ascore=($datos_ascore{$p1}{$nivel_maxima_diferencia}{score})-($datos_ascore{$p2}{$nivel_maxima_diferencia}{score});	
	#return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},$p2."(".$punto1."-".$punto2.")"," ",$ascore,$proporcion);
	my %result;
	my $numero_deshidrataciones=0;
	my $pdh=$p1;
		while($pdh=~s/23//){
			$numero_deshidrataciones++;
		}
	my $numero_fosforilaciones=0;
	my $pph=$p1;
		while($pph=~s/(\(21\)|\(21\,|21\))//){
			$numero_fosforilaciones++;
		}
	$result{'phosphorylations'}=$numero_fosforilaciones;
	$result{'dehydratations'}=$numero_deshidrataciones;
	$result{'peptide'}=$p1;
	$result{'peptide_score'}=$peptide_scores{$p1};
	if($ascore<0.1){
		$ascore=0;
	}
	$result{'ascore'}=[$ascore];
	
	return %result;

}


sub ascore_2_fosforilaciones{
	#print "2 fosfos\n";
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my $masses=$_[2];			#referencia a array de masas
	my $intensities=$_[3];	#referencia a array intensidades
	my $z=$_[4];				#carga
	my $limit_mz_sup=$_[5];			#Limite del maximo valor mz usado (2000?)
	my $limit_mz_inf=$_[6];
	#print "Ascore 2 fosfos\n";
	my @masas=split /\|/,$masses;
	my @intensidades=split /\|/,$intensities;

	
	my @peptidos_posibles=peptide_collection($peptido); 
	#foreach(@peptidos_posibles){
		#print "\t\t",$_,"\n";
	#}
	
	
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z,$limit_mz_sup,$limit_mz_inf)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	

	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}

	#Ordeno peptidos segun Peptide score
	my @peptidos_con_peptides_scores_ordenados=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_peptides_scores_ordenados[0];
	
	#Convierto peptidos en formato [styuv] segun se vaya necesitando
	
	#Detecto puntos de modificacion
	my @p_trozos_peptido_1= split //,conversion_styuv($p1);
	my $f1=0;
	my $f2=0;
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			my $punto=1+$i;
			if($f1){
				$f2=$punto;
			}
			else{
				$f1=$punto;
			}
		}
	}
	
	#Encuentro Peptidos-2 para cada modificacion
	my $peptido2_f1;
	my $peptido2_f2;
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($peptidos_con_peptides_scores_ordenados[$i]);
		if(($p_trozos_peptido_x[(-1+$f1)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)){
			$peptido2_f1=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($peptidos_con_peptides_scores_ordenados[$i]);
		if(($p_trozos_peptido_x[(-1+$f2)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)){
			$peptido2_f2=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}


	#Obtengo niveles maxima diferencia entre Peptide Scores (P1 respecto a peptido2_f1 y peptido2_f2)
	my $nivel_maxima_diferencia1=0;
	my $difss1=0;
	my $nivel_maxima_diferencia2=0;
	my $difss2=0;
	{
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f1}{$i}{score};
			if($dif>$difss1){
				$difss1=$dif;
				$nivel_maxima_diferencia1=$i;
			}
		}
		if($nivel_maxima_diferencia1==0){
			$nivel_maxima_diferencia1=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f2}{$i}{score};
			if($dif>$difss2){
				$difss2=$dif;
				$nivel_maxima_diferencia2=$i;
			}
		}
		if($nivel_maxima_diferencia2==0){
			$nivel_maxima_diferencia2=5;#Esto ocurre cuando los scores del pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
	}
	
	
	#Obtengo fragmentos dentro de los intervalos adecuados
		#1#Encuentro puntos complementarios a las f1 y f2 en cada peptido
		my @p_trozos_peptido_2f1= split //,conversion_styuv($peptido2_f1);
		my @p_trozos_peptido_2f2= split //,conversion_styuv($peptido2_f2);
		my $aa_complement_1;
		for(my $i=0;$i<$#p_trozos_peptido_2f1+1;$i++){
			if($p_trozos_peptido_2f1[$i]=~/[a-z]/){
				unless((1+$i)==$f2){
					$aa_complement_1=1+$i;
				}
			}
		}
		my $aa_complement_2;
		for(my $i=0;$i<$#p_trozos_peptido_2f2+1;$i++){
			if($p_trozos_peptido_2f2[$i]=~/[a-z]/){
				unless((1+$i)==$f1){
					$aa_complement_2=1+$i;
				}
			}
		}
		
		#2#Ordeno los puntos que marcan los intervalos adecuados
			my ($punto11,$punto12,$punto21,$punto22);
			if($f1>$aa_complement_1){
				$punto11=$aa_complement_1;
				$punto12=$f1;
			}
			elsif($f1<$aa_complement_1){
				$punto12=$aa_complement_1;
				$punto11=$f1;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
			
			if($f2>$aa_complement_2){
				$punto21=$aa_complement_2;
				$punto22=$f2;
			}
			elsif($f2<$aa_complement_2){
				$punto22=$aa_complement_2;
				$punto21=$f2;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
		
		
		
		#3#Calculo fragmentos
			#Fragmentos para fosforilacion 1
			my %fragmentos_p1_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p1_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p1_1,$punto11,$punto12)};
			my %fragmentos_p1_2=%{prever_fragmentos($peptido2_f1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p1_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f1),\%fragmentos_p1_2,$punto11,$punto12)};
			#Fragmentos para fosforilacion 2
			my %fragmentos_p2_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p2_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p2_1,$punto21,$punto22)};
			my %fragmentos_p2_2=%{prever_fragmentos($peptido2_f2,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p2_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f2),\%fragmentos_p2_2,$punto21,$punto22)};
		

		##Obtengo datos-ascore: con matches de los fragmentos site-determining
		my $datos_ascore_p1_1=matches_ascore(\%fragmentos_p1_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p1_2=matches_ascore(\%fragmentos_p1_2,\@masas,\@intensidades,$precision,$peptido2_f1);
				
		my $datos_ascore_p2_1=matches_ascore(\%fragmentos_p2_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p2_2=matches_ascore(\%fragmentos_p2_2,\@masas,\@intensidades,$precision,$peptido2_f2);
		
		
		#print Dumper $datos_ascore_p1_1;
		
		
		#print $#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		#exit;
		my $matches1_1=1+$#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my $matches1_2=1+$#{${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my @array1=keys(%fragmentos_p1_1);
		my $size1=1+$#array1;
		my $proporcion1=$matches1_1."//".$size1." (".$matches1_2.")";
		
		my $matches2_1=1+$#{${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my $matches2_2=1+$#{${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my @array2=keys(%fragmentos_p2_1);
		my $size2=1+$#array2;
		my $proporcion2=$matches2_1."//".$size2." (".$matches2_2.")";
		
		
		
		#print "Proporcion1: ",$proporcion1,"\n";
	##	Fosfo1, peptido1
		
		#if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})){
			##print "No definido score p1\n";
			#return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
		#}
		
		#${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{score};
		
		my $ascore1=${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{score}-${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{score};	
		my $ascore2=${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{score}-${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{score};	
		
		#print $p1," -> ";
		#print $ascore1,"\t|\t",$ascore2,"\n";
	my %result;
	my $numero_deshidrataciones=0;
	my $pdh=$p1;
		while($pdh=~s/23//){
			$numero_deshidrataciones++;
		}
	my $numero_fosforilaciones=0;
	my $pph=$p1;
		while($pph=~s/(\(21\)|\(21\,|21\))//){
			$numero_fosforilaciones++;
		}
		
	$result{'phosphorylations'}=$numero_fosforilaciones;
	$result{'dehydratations'}=$numero_deshidrataciones;
	$result{'peptide'}=$p1;
	$result{'peptide_score'}=$peptide_scores{$p1};
	if($ascore1<0.1){
		$ascore1=0;
	}
	if($ascore2<0.1){
		$ascore2=0;
	}
	$result{'ascore'}=[$ascore1,$ascore2];
	
	return %result;	
	#return ($p1,$nivel_maxima_diferencia1." | ".$nivel_maxima_diferencia2,$peptide_scores{$p1},$peptido2_f1."(".$punto11."-".$punto12.")",$peptido2_f2."(".$punto21."-".$punto22.")",$ascore1,$proporcion1,$ascore2,$proporcion2);
		
	
	#my %fragmentos_p2=%{prever_fragmentos($p2,$z)};
	#%fragmentos_p2=%{tratamiento_fragmentos($p2,\%fragmentos_p2,$punto1,$punto2)};
	
	#print $p1,"\n";
	#print "\t$peptido2_f1\n";
	#print "\t$peptido2_f2\n";
	#print "n1: $nivel_maxima_diferencia1\n";
	#print "Puntos 1: $punto11->$punto12\n";
	#print "n2: $nivel_maxima_diferencia2\n";
	#print "Puntos 1: $punto21->$punto22\n\n";	

}


sub ascore_3_fosforilaciones{
	my $peptido=$_[0];			#en formato [styuv]
	my $precision=$_[1];		#en unidades m/z
	my $masses=$_[2];			#referencia a array de masas
	my $intensities=$_[3];	#referencia a array intensidades
	my $z=$_[4];				#carga
	my $limit_mz_sup=$_[5];			
	my $limit_mz_inf=$_[6];
	#print "Ascore 2 fosfos\n";
	my @masas=split /\|/,$masses;
	my @intensidades=split /\|/,$intensities;

	
	my @peptidos_posibles=peptide_collection($peptido); 
	#foreach(@peptidos_posibles){
		#print "\t\t",$_,"\n";
	#}
	
	
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z,$limit_mz_sup,$limit_mz_inf)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	

	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}

	#Ordeno peptidos segun Peptide score
	my @peptidos_con_peptides_scores_ordenados=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_peptides_scores_ordenados[0];
	
	#Convierto peptidos en formato [styuv] segun se vaya necesitando
	
	#Detecto puntos de modificacion
	my @p_trozos_peptido_1= split //,conversion_styuv($p1);
	my $f1=0;
	my $f2=0;
	my $f3=0;
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			my $punto=1+$i;
			if($f2){
				$f3=$punto;
			}
			elsif($f1){
				$f2=$punto;
			}
			else{
				$f1=$punto;
			}
		}
	}


	#Encuentro Peptidos-2 para cada modificacion
	#Los peptidos-2 se escogen con este criterio: se busca un peptido con el mayor peptide-score posible y que
	#	cumpla la condicion que el p1 NO este modificado y que los otros p1(con 3 modificaciones, para el caso de
	# 	buscar para f1, el p1 de f2 y el p1 de f3) SI esten modificados.	
	my $peptido2_f1;
	my $peptido2_f2;
	my $peptido2_f3;
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($peptidos_con_peptides_scores_ordenados[$i]);
		if(($p_trozos_peptido_x[(-1+$f1)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)and($p_trozos_peptido_x[(-1+$f3)]=~/[a-z]/)){
			$peptido2_f1=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($peptidos_con_peptides_scores_ordenados[$i]);
		if(($p_trozos_peptido_x[(-1+$f2)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)and ($p_trozos_peptido_x[(-1+$f3)]=~/[a-z]/)){
			$peptido2_f2=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
	for (my $i=1;$i<$#peptidos_con_peptides_scores_ordenados+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($peptidos_con_peptides_scores_ordenados[$i]);
		if(($p_trozos_peptido_x[(-1+$f3)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)and ($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)){
			$peptido2_f3=$peptidos_con_peptides_scores_ordenados[$i];
			last;
		}
	}
	

	#Obtengo niveles maxima diferencia entre Peptide Scores (P1 respecto a peptido2_f1, peptido2_f2 y peptido2_f3)
	my $nivel_maxima_diferencia1=0;
	my $difss1=0;
	my $nivel_maxima_diferencia2=0;
	my $difss2=0;
	my $nivel_maxima_diferencia3=0;
	my $difss3=0;
	{
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f1}{$i}{score};
			if($dif>$difss1){
				$difss1=$dif;
				$nivel_maxima_diferencia1=$i;
			}
		}
		if($nivel_maxima_diferencia1==0){
			$nivel_maxima_diferencia1=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f2}{$i}{score};
			if($dif>$difss2){
				$difss2=$dif;
				$nivel_maxima_diferencia2=$i;
			}
		}
		if($nivel_maxima_diferencia2==0){
			$nivel_maxima_diferencia2=5;#Esto ocurre cuando los scores del pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f3}{$i}{score};
			if($dif>$difss3){
				$difss3=$dif;
				$nivel_maxima_diferencia1=$i;
			}
		}
		if($nivel_maxima_diferencia3==0){
			$nivel_maxima_diferencia3=5;#Esto ocurre cuando el pep1 y el pep3 son identicos; uso el 5 de forma arbitraria 
		}
	}
	
	
	#Obtengo fragmentos dentro de los intervalos adecuados
		#1#Encuentro puntos complementarios a las f1,f2 y f3 en cada peptido
		#	Por punto complementario se entiende al aminoacido que este modificado en p2 pero NO en p1
		my @p_trozos_peptido_2f1= split //,conversion_styuv($peptido2_f1);
		my @p_trozos_peptido_2f2= split //,conversion_styuv($peptido2_f2);
		my @p_trozos_peptido_2f3= split //,conversion_styuv($peptido2_f3);
		my $aa_complement_1;
		for(my $i=0;$i<$#p_trozos_peptido_2f1+1;$i++){
			if($p_trozos_peptido_2f1[$i]=~/[a-z]/){
				unless(((1+$i)==$f2)or((1+$i)==$f3)){
					$aa_complement_1=1+$i;
				}
			}
		}
		my $aa_complement_2;
		for(my $i=0;$i<$#p_trozos_peptido_2f2+1;$i++){
			if($p_trozos_peptido_2f2[$i]=~/[a-z]/){
				unless(((1+$i)==$f1)or((1+$i)==$f3)){
					$aa_complement_2=1+$i;
				}
			}
		}
		my $aa_complement_3;
		for(my $i=0;$i<$#p_trozos_peptido_2f3+1;$i++){
			if($p_trozos_peptido_2f3[$i]=~/[a-z]/){
				unless(((1+$i)==$f1)or((1+$i)==$f2)){
					$aa_complement_3=1+$i;
				}
			}
		}
		
		
		
		#2#Ordeno los puntos que marcan los intervalos adecuados
		#	Esto significa que ordeno los pares f1-aa_complement1, pares f2-aa_complement2 y 
		#	pares f3-aa_complement3 como pares menor-mayor (sea cual sea el de p1 y p2)
		   
			my ($punto11,$punto12,$punto21,$punto22,$punto31,$punto32);
			if($f1>$aa_complement_1){
				$punto11=$aa_complement_1;
				$punto12=$f1;
			}
			elsif($f1<$aa_complement_1){
				$punto12=$aa_complement_1;
				$punto11=$f1;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
			
			if($f2>$aa_complement_2){
				$punto21=$aa_complement_2;
				$punto22=$f2;
			}
			elsif($f2<$aa_complement_2){
				$punto22=$aa_complement_2;
				$punto21=$f2;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
			if($f3>$aa_complement_3){
				$punto31=$aa_complement_3;
				$punto32=$f3;
			}
			elsif($f3<$aa_complement_3){
				$punto32=$aa_complement_3;
				$punto31=$f3;
			}
			else{
				print "errorrrrrrrrrrrrrrr!!!!\n";
			}
		
		#print conversion_styuv($p1),"\n\n";
		#print "pep2_f1: ",conversion_styuv($peptido2_f1),"\n";
		#print "max diff 1: ",$nivel_maxima_diferencia1,"\n";
		#print "punto complementario f1: ",$aa_complement_1,"\n";
		#print "punto 11: $punto11 , punto 12: $punto12 \n";
		#print "pep2_f2: ",conversion_styuv($peptido2_f2),"\n";
		#print "max diff 2: ",$nivel_maxima_diferencia2,"\n";
		#print "punto complementario f2: ",$aa_complement_2,"\n";
		#print "punto 21: $punto21 , punto 22: $punto22 \n";
		#print "pep2_f3: ",conversion_styuv($peptido2_f3),"\n";
		#print "max diff 3: ",$nivel_maxima_diferencia3,"\n";
		#print "punto complementario f3: ",$aa_complement_3,"\n";
		#print "punto 31: $punto31 , punto 32: $punto32 \n";
		#print "------------------------------------------------\n\n";
		#return;
		
		#3#Calculo fragmentos
			#Fragmentos para fosforilacion 1
			my %fragmentos_p1_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p1_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p1_1,$punto11,$punto12)};
			my %fragmentos_p1_2=%{prever_fragmentos($peptido2_f1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p1_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f1),\%fragmentos_p1_2,$punto11,$punto12)};
			#Fragmentos para fosforilacion 2
			my %fragmentos_p2_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p2_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p2_1,$punto21,$punto22)};
			my %fragmentos_p2_2=%{prever_fragmentos($peptido2_f2,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p2_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f2),\%fragmentos_p2_2,$punto21,$punto22)};
			#Fragmentos para fosforilacion 3
			my %fragmentos_p3_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p3_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p3_1,$punto31,$punto32)};
			my %fragmentos_p3_2=%{prever_fragmentos($peptido2_f3,$z,$limit_mz_sup,$limit_mz_inf)};
			%fragmentos_p3_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f3),\%fragmentos_p3_2,$punto31,$punto32)};
			
			
			

		##Obtengo datos-ascore: con matches de los fragmentos site-determining
		my $datos_ascore_p1_1=matches_ascore(\%fragmentos_p1_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p1_2=matches_ascore(\%fragmentos_p1_2,\@masas,\@intensidades,$precision,$peptido2_f1);
				
		my $datos_ascore_p2_1=matches_ascore(\%fragmentos_p2_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p2_2=matches_ascore(\%fragmentos_p2_2,\@masas,\@intensidades,$precision,$peptido2_f2);
		
		my $datos_ascore_p3_1=matches_ascore(\%fragmentos_p3_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p3_2=matches_ascore(\%fragmentos_p3_2,\@masas,\@intensidades,$precision,$peptido2_f3);
		
		#print Dumper $datos_ascore_p1_1;
		
		
		#print $#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		#exit;
		my $matches1_1=1+$#{${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my $matches1_2=1+$#{${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{tags_match_intervalo}};
		my @array1=keys(%fragmentos_p1_1);
		my $size1=1+$#array1;
		my $proporcion1=$matches1_1."//".$size1." (".$matches1_2.")";
		
		my $matches2_1=1+$#{${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my $matches2_2=1+$#{${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{tags_match_intervalo}};
		my @array2=keys(%fragmentos_p2_1);
		my $size2=1+$#array2;
		my $proporcion2=$matches2_1."//".$size2." (".$matches2_2.")";
		
		my $matches3_1=1+$#{${$datos_ascore_p3_1}{$nivel_maxima_diferencia3}{tags_match_intervalo}};
		my $matches3_2=1+$#{${$datos_ascore_p3_2}{$nivel_maxima_diferencia3}{tags_match_intervalo}};
		my @array3=keys(%fragmentos_p3_1);
		my $size3=1+$#array3;
		my $proporcion3=$matches3_1."//".$size3." (".$matches3_2.")";
		
		#print "Proporcion1: ",$proporcion1,"\n";
	##	Fosfo1, peptido1
		
		#if((not defined $datos_ascore{$p1}{$nivel_maxima_diferencia}{score})){
			##print "No definido score p1\n";
			#return ($p1,$nivel_maxima_diferencia,$peptide_scores{$p1},0,$proporcion);
		#}
		
		#${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{score};
		
		my $ascore1=${$datos_ascore_p1_1}{$nivel_maxima_diferencia1}{score}-${$datos_ascore_p1_2}{$nivel_maxima_diferencia1}{score};	
		my $ascore2=${$datos_ascore_p2_1}{$nivel_maxima_diferencia2}{score}-${$datos_ascore_p2_2}{$nivel_maxima_diferencia2}{score};	
		my $ascore3=${$datos_ascore_p3_1}{$nivel_maxima_diferencia3}{score}-${$datos_ascore_p3_2}{$nivel_maxima_diferencia3}{score};
		#print $p1," -> ";
		#print $ascore1,"\t|\t",$ascore2,"\n";
		
		my %result;
		my $numero_deshidrataciones=0;
		my $pdh=$p1;
			while($pdh=~s/23//){
				$numero_deshidrataciones++;
			}
		my $numero_fosforilaciones=0;
		my $pph=$p1;
			while($pph=~s/(\(21\)|\(21\,|21\))//){
				$numero_fosforilaciones++;
			}
			
		$result{'phosphorylations'}=$numero_fosforilaciones;
		$result{'dehydratations'}=$numero_deshidrataciones;
		$result{'peptide'}=$p1;
		$result{'peptide_score'}=$peptide_scores{$p1};
		if($ascore1<0.1){
			$ascore1=0;
		}
		if($ascore2<0.1){
			$ascore2=0;
		}
		if($ascore3<0.1){
			$ascore3=0;
		}
		$result{'ascore'}=[$ascore1,$ascore2,$ascore3];
		
		return %result;	
		
		
		#return ($p1,$nivel_maxima_diferencia1." | ".$nivel_maxima_diferencia2." | ".$nivel_maxima_diferencia3,$peptide_scores{$p1},$peptido2_f1."(".$punto11."-".$punto12.")",$peptido2_f2."(".$punto21."-".$punto22.")",$peptido2_f3."(".$punto31."-".$punto32.")",$ascore1,$proporcion1,$ascore2,$proporcion2,$ascore3,$proporcion3);
		
	
	#my %fragmentos_p2=%{prever_fragmentos($p2,$z)};
	#%fragmentos_p2=%{tratamiento_fragmentos($p2,\%fragmentos_p2,$punto1,$punto2)};
	
	#print $p1,"\n";
	#print "\t$peptido2_f1\n";
	#print "\t$peptido2_f2\n";
	#print "n1: $nivel_maxima_diferencia1\n";
	#print "Puntos 1: $punto11->$punto12\n";
	#print "n2: $nivel_maxima_diferencia2\n";
	#print "Puntos 1: $punto21->$punto22\n\n";	

}

sub peptide_collection{
	my $pep=$_[0];
	#Asumo dos cosas:
	#	1/Solo hay un itraq por cada aminoacido
	#	2/Los itraq no cambian de sitio
	my @posibles_puntos_fosforilacion;
	my $numero_fosforilaciones=0;
	my @posibles_puntos_DH;
	my $numero_DH=0;
	my $numero_itraq=0;

	my @aas;
	my @mods;
	my $aa_num=0;

	while($pep=~/\w/){
		my $aa;
		my $mod;
		if($pep=~s/^(\w)\((\w+)\)//){
			$aa=$1;
			$mod=$2;
			push @aas,$aa;
			push @mods,$mod;
			#Calculo el numero de fosfos/itraqs/DH
			if($mod==21){
				$numero_fosforilaciones++;
			}
			elsif($mod==23){
				$numero_DH++;
			}
			elsif($mod==737){
				$numero_itraq++;
				if($aa_num>0){#Si hay itraq en un aa no N-terminal, no se puede modificar con nada mas
					$aa_num++;
					next;
				}
			}
			#Recopilo los posibles puntos de fosfo/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		else{
			$pep=~s/^(\w)//;
			$aa=$1;
			$mod='';
			push @aas,$aa;
			push @mods,$mod;
			#Recopilo los posibles puntos de fosfo/itraq/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		$aa_num++;
	}
	#print "posibles puntos DH: ";
	#foreach(@posibles_puntos_DH){
		#print $_,",";
	#}
	#print "\n";
		
	
	#Obtengo las diferentes combinaciones de fosforilados, deshidratados y itraq
	my @combinaciones_fosfo;
	my @combinaciones_DH;
	my @combinaciones_itraq;
	
	
	my $combinat1 = Math::Combinatorics->new(count => $numero_fosforilaciones,data => [@posibles_puntos_fosforilacion],);
	while(my @combo1 = $combinat1->next_combination){
		push @combinaciones_fosfo,\@combo1;
	}
	my $combinat2 = Math::Combinatorics->new(count => $numero_DH,data => [@posibles_puntos_DH],);
		while(my @combo2 = $combinat2->next_combination){
			push @combinaciones_DH,\@combo2;
	}

	#print Dumper @combinaciones_fosfo,"\n";
	#print "\n\n\n";
	#print Dumper @combinaciones_DH,"\n";
	
	my @combinaciones_peptidos;
	if(@combinaciones_fosfo and @combinaciones_DH){
		#print "Hay phospho, DH y itraq\t";
		foreach my $fos(@combinaciones_fosfo){
				my %usados;
				my @fos=@$fos;
				foreach my $f(@fos){
					$usados{$f}=1;
				}
				LABEL1:
				foreach my $dh(@combinaciones_DH){
					my @dhs=@$dh;
					foreach my $d(@dhs){
						if(defined $usados{$d}){
							next LABEL1;
						}
					}
					#Genero un nuevo array de modificaciones para construir el peptido modificado
					my @mods_new;
					for(my $i=0;$i<$#mods+1;$i++){#Recupero las modificaciones no phospho o DH (itraq,MetOx,...)
						if($mods[$i]){
							unless(($mods[$i]==21)or($mods[$i]==23)){
								$mods_new[$i]='('.$mods[$i].')';
							}
						}
					}
					foreach(@fos){
						if($mods_new[$_]){
							my $t=$mods_new[$_];
							$t=~s/[\(,\)]//g;
							$mods_new[$_]='('.$t.',21)';
							
						}
						else{
							$mods_new[$_]='(21)';
						}
					}
					foreach(@dhs){
						if($mods_new[$_]){
							my $t=$mods_new[$_];
							$t=~s/[\(,\)]//g;
							$mods_new[$_]='('.$t.',23)';
							
						}
						else{
							$mods_new[$_]='(23)';
						}
					}
					#Genero el peptido
					my $peptide_string='';
					for(my $i=0;$i<$#aas+1;$i++){
						if($mods_new[$i]){
							$peptide_string=$peptide_string.$aas[$i].$mods_new[$i];
						}
						else{
							$peptide_string=$peptide_string.$aas[$i];
						}
					}
					push @combinaciones_peptidos,$peptide_string;
				}
		}
	}
	elsif(@combinaciones_fosfo){
		#print "solo fosfo\t";
		foreach my $fos(@combinaciones_fosfo){
				my @fos=@$fos;
				#Genero un nuevo array de modificaciones para construir el peptido modificado
					my @mods_new;
					for(my $i=0;$i<$#mods+1;$i++){#Recupero las modificaciones no phospho o DH (itraq,MetOx,...)
						if($mods[$i]){
							unless(($mods[$i]==21)or($mods[$i]==23)){
								$mods_new[$i]='('.$mods[$i].')';
							}
						}
					}
					foreach(@fos){
						if($mods_new[$_]){
							my $t=$mods_new[$_];
							$t=~s/[\(,\)]//g;
							$mods_new[$_]='('.$t.',21)';
							
						}
						else{
							$mods_new[$_]='(21)';
						}
					}
					#Genero el peptido
					my $peptide_string='';
					for(my $i=0;$i<$#aas+1;$i++){
						if($mods_new[$i]){
							$peptide_string=$peptide_string.$aas[$i].$mods_new[$i];
						}
						else{
							$peptide_string=$peptide_string.$aas[$i];
						}
					}
					push @combinaciones_peptidos,$peptide_string;
		
		}
	}
	elsif(@combinaciones_DH){
		#print "solo DH\t";
		foreach my $dhs(@combinaciones_DH){
				my @dhs=@$dhs;
				#Genero un nuevo array de modificaciones para construir el peptido modificado
					my @mods_new;
					for(my $i=0;$i<$#mods+1;$i++){#Recupero las modificaciones no phospho o DH (itraq,MetOx,...)
						if($mods[$i]){
							unless(($mods[$i]==21)or($mods[$i]==23)){
								$mods_new[$i]='('.$mods[$i].')';
							}
						}
					}
					foreach(@dhs){
						if($mods_new[$_]){
							my $t=$mods_new[$_];
							$t=~s/[\(,\)]//g;
							$mods_new[$_]='('.$t.',23)';
							
						}
						else{
							$mods_new[$_]='(23)';
						}
					}
					#Genero el peptido
					my $peptide_string='';
					for(my $i=0;$i<$#aas+1;$i++){
						if($mods_new[$i]){
							$peptide_string=$peptide_string.$aas[$i].$mods_new[$i];
						}
						else{
							$peptide_string=$peptide_string.$aas[$i];
						}
					}
					push @combinaciones_peptidos,$peptide_string;
		}
		
	}
	else{
		push @combinaciones_peptidos,$_[0];
	}

	return @combinaciones_peptidos;	
}

sub prever_fragmentos{
	# Module to predict a spectra using a sequence. It uses the 20 aa and allows phosphorilations on STY (+80:s,t,y),
	#water loss on ST (-18:u,v), and methionine oxidation (+18:m). It produces a hash containing the {label}->mass.
	#It calculates the y and b series for the charge state, 
	#ARGUMENTS: there are 3 args:
	#       - $_[0]-> sequence
	#       - $_[1]-> charge state (1,2 or 3)
	#       - $_[2]-> mz upper limit for fragment generation
	#RETURN: a hash with {label}->mass.
	my $peptide=$_[0];
	my $z=$_[1];
	my $limit_mz_sup=$_[2];
	my $limit_mz_inf=$_[3];
	my %masas_aa = (
	A=>71.03711,V=>99.06841,L=>113.08406,I=>113.08406,P=>97.05276,M	=>131.0405,
	F=>147.0684,W=>186.07932,G=>57.02146,S=>87.03203,T=>101.04768,C=>160.03069,
	Y=>163.06332,N=>114.04293,Q=>128.05858,D=>115.02694,E=>129.04259,H=>137.05891,
	K=>128.09496,R=>156.1011);
	#la C esta CAM!, la considero PTM fija
	
	my @aas;
	my @mods;
	my @masas_aas;

	while($peptide=~/\w/){
		my $aa;
		my $mod;
		if($peptide=~s/^(\w)\((\w+\,\w+)\)//){
			$aa=$1;
			$mod=$2;
			push @aas,$aa;
			push @mods,$mod;
		}
		elsif($peptide=~s/^(\w)\((\w+)\)//){
			$aa=$1;
			$mod=$2;
			push @aas,$aa;
			push @mods,$mod;
		}
		else{
			$peptide=~s/^(\w)//;
			$aa=$1;
			$mod='';
			push @aas,$aa;
			push @mods,$mod;
		}
	}
	my %masas_mod=(21=>79.9799,23=>-18.010565,737=>229.162932,35=>15.994915);
	for(my $i=0;$i<$#aas+1;$i++){
		if($mods[$i]){
			if($mods[$i]=~/\,/){
				my($m1,$m2)=split /\,/,$mods[$i];
				$masas_aas[$i]=$masas_mod{$m1}+$masas_mod{$m2}+$masas_aa{$aas[$i]};
			}
			else{
				$masas_aas[$i]=$masas_mod{$mods[$i]}+$masas_aa{$aas[$i]};
				#print $mods[$i],"\n";
			}
		}
		else{
			$masas_aas[$i]=$masas_aa{$aas[$i]};
		}
	}
	#foreach(@masas_aas){
		#print $_,"\n";
	#}
	#print "\n\n";
	##exit;
	
	#La carga maxima permitida es 1 para z=1 o 2, 
	my $carga_maxima_permitida_fragmentos;
	if($z==1){
		$carga_maxima_permitida_fragmentos=1;
	}
	elsif($z==2){
		$carga_maxima_permitida_fragmentos=1;
	}
	else{
		$carga_maxima_permitida_fragmentos=$z-1;
	}


	my $acumulado=0;
	my (@serie_y_1,@serie_b_1,@serie_y_2,@serie_b_2,@serie_y_3,@serie_b_3);

	for (my $i=0;$i<$#aas+1;$i++){
		$serie_b_1[$i]=$masas_aas[$i]+$acumulado+1;
		$acumulado=$acumulado+$masas_aas[$i];
		}
	my $M=$acumulado+19;
	for(my $i=0;$i<$#aas;$i++){
		$serie_y_1[$i]=$M-$serie_b_1[($#aas-$i)]+1;
		}
	$serie_y_1[$#aas]=$M-$serie_b_1[0]+1;
	if($_[1]>1){
		for (my $i=0;$i<$#aas+1;$i++){
			$serie_y_2[$i]=($serie_y_1[$i]+1)/2;
			if (defined $serie_b_1[$i]){
				$serie_b_2[$i]=($serie_b_1[$i]+1)/2;
				}
			}
		if($_[1]>2){
			for (my $i=0;$i<$#aas+1;$i++){
				$serie_y_3[$i]=($serie_y_1[$i]+1)/3;
				$serie_b_3[$i]=($serie_b_1[$i]+1)/3;
			}		
		}	    	
}


	my %fragmentos;

	##Prueba de fragmentos especiales
	#my $M_perdida98=($M+$_[1]-98)/$_[1];
	#$fragmentos{"esp1"}=sprintf("%.3f",$M_perdida98);	
	#$fragmentos{"esp1b"}=sprintf("%.3f",$M_perdida98/$_[1]);	
	#my $M_perdidaagua=($M+$_[1]-18)/$_[1];
	#$fragmentos{"esp2"}=sprintf("%.3f",$M_perdidaagua);
	#$fragmentos{"esp2b"}=sprintf("%.3f",$M_perdidaagua/$_[1]);


	for (my $i=0;$i<$#aas+1;$i++){
		my $j=$i+1;
		if ($i<($#aas)){											# El b(n) bo existe
			my $num=$j;
			if($num<10){$num='0'.$num};
			$fragmentos{"b$num+1"}=sprintf("%.3f",$serie_b_1[$i]);
			if(($_[1]>1)and($carga_maxima_permitida_fragmentos>1)){													# Si carga>1
				$fragmentos{"b$num+2"}=sprintf("%.3f",$serie_b_2[$i]);
				if(($_[1]>2)and($carga_maxima_permitida_fragmentos>2)){												# Si carga>2
					$fragmentos{"b$num+3"}=sprintf("%.3f",$serie_b_3[$i]);
				}
			}	
		}
		if(($i>0)){	
			my $num=$i;
			if($num<10){$num='0'.$num};														# El y(0) bo existe
			$fragmentos{"y$num+1"}=sprintf("%.3f",$serie_y_1[$i]);
			 if(($_[1]>1)and($carga_maxima_permitida_fragmentos>1)){													# Si carga>1
				$fragmentos{"y$num+2"}=sprintf("%.3f",$serie_y_2[$i]);
				if(($_[1]>2)and($carga_maxima_permitida_fragmentos>2)){												# Si carga>2
					$fragmentos{"y$num+3"}=sprintf("%.3f",$serie_y_3[$i]);	
				}	
			 }
		}	
	}
	#Solo se consideran los fragmentos entre limit_mz_sup e _inf
	foreach my $frag(keys %fragmentos){
		if(($fragmentos{$frag}>$limit_mz_sup)or($fragmentos{$frag}<$limit_mz_inf)){
			delete $fragmentos{$frag};
		}
	}

	return \%fragmentos;
}


sub matches_10niveles_series_yb_completas{
	my %fragmentos=%{$_[0]};
	my @masas=@{$_[1]};
	my @intensidades=@{$_[2]};
	my $fragment_error=$_[3];
	my $peptido=$_[4];
	
	my $numero_fragmentos=0;
	foreach(keys %fragmentos){
		$numero_fragmentos++;
		#print $_,"\t",$fragmentos{$_},"\n";
	}
	
	
	my $num_picos_cada_100;
	my %datos_numeros_picos;
	my @rangos_masas;
	my ($rango_masas_inf,$rango_masas_sup)=(200,2000);
	for(my $i=$rango_masas_inf;$i<$rango_masas_sup;$i=$i+100){
		push @rangos_masas,$i;
	}
	for($num_picos_cada_100=1;$num_picos_cada_100<11;$num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias=0;
		my $num_picos_dentro_intervalo=0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior=$intervalo-100;
			for(my $i=0;$i<$#masas+1;$i++){
				if(($masas[$i]>$intervalo_inferior)and($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]}=$intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			if($num_picos_dentro_intervalo==0){#ningun pico experimental en el intervalo
										       #No hay picos en el intervalo\n";
				next;
			}
			elsif($num_picos_dentro_intervalo==1){ #un solo pico experimental en el intervalo
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			else{#mas de 1 pico experimental en el intervalo
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
				
			}
			
			#	Aqui se realiza el matching entre masas teoricas y experimentales
			my $masas_testadas=0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}
				#$numero_fragmentos=> trials
				#$correspondencias=>  successes
				#$prob=> numero picos cada 100 / 100

		my $probablity=($num_picos_cada_100/100);
		my $binomial_acum_rango=binomial($numero_fragmentos,$correspondencias,$probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score=-10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score=0;
		}
		
		$datos_numeros_picos{$num_picos_cada_100}{masas_exp_match_intervalo}=\@masas_exp_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{masas_teor_match_intervalo}=\@masas_teor_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{tags_match_intervalo}=\@tags_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{score}=$score;
	}
	
	return \%datos_numeros_picos;
}

sub binomial{
	#Binomial distribution
	#($n,$k,$p)->$k successes in $n tries, given a probability of $p	
	my $n=$_[0];
	my $k=$_[1];
	my $p=$_[2];
	my $binomial_acumulada=0;
	for(my $i=$k; $i<$n+1;$i++){
		my $binomial=Math::NumberCruncher::Binomial($n,$i,$p);
		$binomial_acumulada= $binomial_acumulada+$binomial;
	}
	return $binomial_acumulada;
}

sub conversion_styuv{
	my $pep=$_[0];
	#cambio fosfo o DH simples
	$pep=~s/S\(21\)/s/g;
	$pep=~s/T\(21\)/t/g;
	$pep=~s/Y\(21\)/y/g;
	$pep=~s/S\(23\)/u/g;
	$pep=~s/T\(23\)/v/g;
	#cambio fosfo o DH combinados con iTraq en Nt
	$pep=~s/S\(737\,21\)/s/;
	$pep=~s/T\(737\,21\)/t/;
	$pep=~s/Y\(737\,21\)/y/;
	$pep=~s/S\(21\,737\)/s/;
	$pep=~s/T\(21\,737\)/t/;
	$pep=~s/Y\(21\,737\)/y/;
	$pep=~s/S\(737\,23\)/u/;
	$pep=~s/T\(737\,23\)/v/;
	$pep=~s/S\(23\,737\)/u/;
	$pep=~s/T\(23\,737\)/v/;
	#Elimino todas las demas modificaciones
	$pep=~s/\(\d+\)//g;
	return $pep;
}

sub tratamiento_fragmentos{
	#argumentos-> $_[0]:%fragmentos,$_[1]:$peptido,$_[2]:primer aa intervalo, $_[3]: ultimo aa intervalo;
	########		Tratamiento fragmentos		######
	#Aqui genero un nuevo hash para poner los fragmentos que estan entre el primero y ultimo aa que se pueda fosforilar
		
		my $peptido=$_[0];
		my %fragmentos=%{$_[1]};
		#print Dumper %fragmentos;
		#exit;
		my $primer_aa_fosforilable=$_[2];
		my $ultimo_aa_fosforilable=$_[3];
		my %fragmentos_filtrados;
	
		my $largo_pep= length $peptido;
		
		##eliminacion de fragmentos especiales
		#foreach my $frag (keys %fragmentos){
			#if($frag=~/^esp/){
				#$fragmentos_filtrados{$frag}=$fragmentos{$frag};
				#delete $fragmentos{$frag};
			#}
		#}
		
		
		foreach my $frag (keys %fragmentos){
			my $frag2=(split /\+/,$frag)[0];
			if($frag2=~m/y/){
				$frag2=~s/y//;		#para las y, valen los fragmentos entre [largo-ultimo+1 , largo-primero]
				if(($frag2>=($largo_pep-$ultimo_aa_fosforilable+1))and($frag2<=($largo_pep-$primer_aa_fosforilable))){
					$fragmentos_filtrados{$frag}=$fragmentos{$frag};
				}
			}
			elsif($frag2=~m/b/){	#para las b, valen los fragmentos entre [primero , ultimo-1]
				$frag2=~s/b//;
				if(($frag2>=($primer_aa_fosforilable))and($frag2<=($ultimo_aa_fosforilable-1))){
					$fragmentos_filtrados{$frag}=$fragmentos{$frag};
				}
			}
		}
		
		
		%fragmentos=();
		%fragmentos = %fragmentos_filtrados;#Uso solo fragmentos en zona de interes		
	#######		fin tratamiento fragmentos
	return \%fragmentos;
}


sub matches_ascore{
	my %fragmentos=%{$_[0]};
	my @masas=@{$_[1]};
	my @intensidades=@{$_[2]};
	my $fragment_error=$_[3];
	my $peptido=$_[4];
	
	
	my $numero_fragmentos=0;
	foreach(keys %fragmentos){
		$numero_fragmentos++;
		#print $_,"\t",$fragmentos{$_},"\n";
	}
	
	
	my $num_picos_cada_100;
	my %datos_numeros_picos;
	my @rangos_masas;
	my ($rango_masas_inf,$rango_masas_sup)=(200,2000);
	for(my $i=$rango_masas_inf;$i<$rango_masas_sup;$i=$i+100){
		push @rangos_masas,$i;
	}
	for($num_picos_cada_100=1;$num_picos_cada_100<11;$num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias=0;
		my $num_picos_dentro_intervalo=0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior=$intervalo-100;
			for(my $i=0;$i<$#masas+1;$i++){
				if(($masas[$i]>$intervalo_inferior)and($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]}=$intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			if($num_picos_dentro_intervalo==0){#ningun pico experimental en el intervalo
										       #No hay picos en el intervalo\n";
				next;
				#####????????????????????????
			}
			elsif($num_picos_dentro_intervalo==1){ #un solo pico experimental en el intervalo
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			else{#mas de 1 pico experimental en el intervalo
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
				#my $kk=1+$#masas_ordenadas_por_intensidad;
				#if($kk>0){
					#print "num picos cada 100: $num_picos_cada_100 , Intervalo: $intervalo , picos: ".$kk."\n";
					#foreach(@masas_ordenadas_por_intensidad){print $_," | "}
					#print "\n";				
				#}
				
			}
			
			#	Aqui se realiza el matching entre masas teoricas y experimentales
			my $masas_testadas=0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}
				#$numero_fragmentos=> trials
				#$correspondencias=>  successes
				#$prob=> numero picos cada 100 / 100

		my $probablity=($num_picos_cada_100/100);
		my $binomial_acum_rango=binomial($numero_fragmentos,$correspondencias,$probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score=-10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score=0;
		}
		
		$datos_numeros_picos{$num_picos_cada_100}{masas_exp_match_intervalo}=\@masas_exp_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{masas_teor_match_intervalo}=\@masas_teor_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{tags_match_intervalo}=\@tags_match_intervalo;
		$datos_numeros_picos{$num_picos_cada_100}{score}=$score;
	}
	
	return \%datos_numeros_picos;



}


